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Abstract 

An  algorithm  is  introduced  for  the  rapid  evaluation  at  appropriately  chosen  nodes 
on  the  two-dimensional  sphere  S 2  in  R3  of  functions  specified  by  their  spherical  har¬ 
monic  expansions  (known  as  the  inverse  spherical  harmonic  transform),  and  for  the 
evaluation  of  the  coefficients  in  spherical  harmonic  expansions  of  functions  specified 
by  their  values  at  appropriately  chosen  points  on  S 2  (known  as  the  forward  spherical 
harmonic  transform).  The  procedure  is  numerically  stable  and  requires  an  amount  of 
CPU  time  proportional  to  iV(log  N)  log(l/e),  where  N  is  the  number  of  nodes  in  the 
discretization  of  S'2,  and  e  is  the  precision  of  computations.  The  performance  of  the 
algorithm  is  illustrated  via  several  numerical  examples. 


1  Introduction 

Spherical  harmonic  expansions  are  a  well-understood  and  widely  used  tool  of  applied  math¬ 
ematics;  they  are  encountered,  inter  alia,  in  weather  and  climate  modeling,  in  the  repre¬ 
sentation  of  gravitational,  topographic,  and  magnetic  data  in  geophysics,  in  the  numerical 
solution  of  certain  partial  differential  equations,  etc.  The  role  of  spherical  harmonic  expan¬ 
sions  in  the  solution  of  the  Laplace  equation  in  three  dimensions  is  similar  to  the  role  played 
by  Fourier  series  expansions  in  two  dimensions. 

The  spherical  harmonic  expansion  of  a  function  /  in  £2(S2)  is  the  series  of  the  form 

OO  l 

m  fiTE  «r  hm|(  cose)  (i) 

1=0  m=—l 

where  ( 6 ,  y)  are  the  standard  spherical  coordinates  on  the  two-dimensional  sphere  S 2  in  R3, 
0  <  9  <  Ti  and  0  <  ip  <  2i r,  and  P/n  is  the  associated  Legendre  function  of  degree  l  and  order 
m.  While  the  functions  {P\'n\cos6)  eimd}  constitute  a  basis  of  £2(S2)  that  is  orthogonal, 
that  is, 

£iPt‘'(x)Plm\x)dx  =  0  (2) 
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when  l  ^  k ,  their  norms  are  not  equal  to  1;  in  fact,  they  are  so  badly  normalized  as  to  be 
virtually  unusable  in  numerical  calculations  (see  Subsection  2.1  below  for  a  detailed  discus¬ 
sion  of  the  associated  Legendre  functions).  Therefore,  it  is  customary  to  replace  expansions 
of  the  form  (1)  with  expansions  of  the  form 

CO  l 

m<p)  =  £  £  «r  P,  '  (cos  9)  (3) 

1=0  m=—l 

where  p\‘n '  denotes  the  normalized  version  of  the  associated  Legendre  function  P\m\  dehned 
on  [—1,  1]  via  the  formula 


p\m\(x) 


21  +  1  (l  —  \m\)\ 
2  (l  +  \m\)\ 


plm+), 


(4) 


so  that 

J^(p'r'(x)f  dx=  1.  (5) 

In  numerical  practice,  the  series  (3)  is  truncated  after  a  finite  number  of  terms,  leading 
to  expressions  of  the  form 

N  1 

m  E  «r d  (cos 9) e4’"*’;  (6) 

/— 0  m=—l 

(6)  is  viewed  as  an  approximation  to  the  function  /,  and  N  is  called  the  order  of  the 
expansion  (6).  Obviously,  the  expansion  (6)  contains  (N  +  l)2  terms;  the  order  N  required 
to  obtain  a  prescribed  accuracy  of  the  approximation  is  determined  by  the  complexity  of  the 
function  /. 

Frequently,  the  need  arises  to  evaluate  the  coefficients  in  an  expansion  of  the  form  (6) 
for  a  function  /  given  by  a  table  of  its  values  at  a  collection  of  appropriately  chosen  nodes 
on  S 2;  conversely,  given  the  coefficients  in  (6),  one  often  needs  to  evaluate  /  at  a  collection 
of  points  on  S2.  The  former  is  usually  called  the  forward  spherical  harmonic  transform,  and 
the  latter  is  known  as  the  inverse  spherical  harmonic  transform.  A  standard  discretization 
of  S2  is  the  “tensor  product,”  consisting  of  all  pairs  of  the  form  (dfc,  <■ Pj ),  with  equispaeed 
nodes  do,  ffi,  . . . ,  dv-i,  9 ^  discretizing  the  interval  [0,  7r],  dehned  by  the  formula 


n(k  +  1/2) 
h=  N  + 1 


(7) 


and  equispaeed  nodes  +0,  +i,  . . . ,  +2N-1,  +2N  discretizing  the  interval  [0,  27t],  dehned  by  the 
formula 


+ 1/2) 

2N  +  1 


(8) 


This  leads  immediately  to  numerical  schemes  for  both  the  forward  and  inverse  spherical 
harmonic  transforms  costing  0(N3)  operations.  Indeed,  given  a  function  /  dehned  on  S2  by 
the  formula  (6),  one  can  rewrite  (6)  in  the  form 


/(0,  <P) 


N  N  i  i 

£  eim *  Y,  aTPi  (cos 9). 

m.=—N  l=\m\ 


(9) 
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For  a  fixed  value  of  9 ,  each  of  the  inner  sums  in  (9)  contains  no  more  than  N  +  1  terms, 
and  there  are  2N+1  such  sums  (one  for  each  value  of  m);  since  the  inverse  spherical  harmonic 
transform  involves  N  +  1  values  6q,  9\ ,  . . . ,  9n- i,  #/v,  the  cost  of  evaluating  all  inner  sums 
in  (9)  is  0(N3).  Once  all  inner  sums  have  been  evaluated,  evaluation  of  each  outer  sum 
costs  O(N)  operations  (since  each  of  them  contains  2N  +  1  terms),  and  there  are  0(N2) 
such  sums  to  be  evaluated,  leading  to  0(N3)  CPU  time  requirements  for  the  evaluation  of 
all  outer  sums  in  (9).  The  cost  of  the  evaluation  of  the  whole  inverse  spherical  harmonic 
transform  (in  the  form  (9))  is  the  sum  of  the  costs  for  the  inner  and  outer  sums,  and  is 
also  0(N3)]  a  virtually  identical  calculation  shows  that  the  cost  of  evaluating  of  the  forward 
harmonic  transform  is  also  0(N3). 

A  trivial  modification  of  the  scheme  described  in  the  preceding  paragraph  uses  the  Fast 
Fourier  Transform  to  evaluate  the  outer  sums  in  (9),  roughly  halving  the  CPU  time  require¬ 
ments  of  the  whole  procedure.  Several  other  considerations  (see,  for  example,  [2])  can  be  used 
to  reduce  the  CPU  time  requirements  by  a  further  factor  of  four  or  so,  but  there  is  no  simple 
trick  for  reducing  the  asymptotic  CPU  time  requirements  of  the  whole  spherical  harmonic 
transform  (either  forward  or  inverse)  below  N3.  In  this  paper,  we  introduce  algorithms 
for  both  forward  and  inverse  spherical  harmonic  transforms  with  CPU  time  requirements 
proportional  to  N2(logN)  log(l/e),  where  £  is  the  precision  of  computations. 

The  algorithm  of  this  paper  is  a  procedure  for  the  rapid  evaluation  of  the  inner  sums  in 
expressions  of  the  form  (9).  It  is  based  principally  on  two  observations,  as  follows: 

1.  The  differential  equations  defining  the  functions  Pjn  with  arbitrary  positive  integer  m 

— 1  — 2 

are  very  close  to  the  differential  equations  defining  the  functions  Pi  and  Pl . 

2.  There  exist  fast  algorithms  for  decomposing  functions  into  and  reconstructing  functions 
from  sums  of  the  forms 

N 

f(x)  =  Y,PiP1i(x)>  (10) 

1=0 

N 

f(x)  =  YJp2iP2i{x)-  (ii) 

1=0 

We  use  the  connections  between  the  functions  Pj'  with  arbitrary  positive  integer  m 
— 1  — 2 

and  the  functions  Pl  and  P{ ,  to  apply  rapidly  to  arbitrary  vectors  the  matrices  converting 
between  expansions  of  the  forms  (10)  and  (11)  and  expansions  of  the  form 

N 

/W  =  yPrAW-  (i2) 

1=0 

This  step  utilizes  the  observation  made  in  [4],  that  the  N  x  N  matrix  of  eigenvectors  of 
the  sum  of  a  diagonal  matrix  and  a  semiseparable  matrix  (see  Subsection  2.4  below  for  the 
definition  of  a  semiseparable  matrix)  can  be  applied  to  an  arbitrary  vector  of  length  N  for  a 
cost  proportional  to  N(logN)  log(l / £:)  operations,  where  £  is  the  precision  of  computations. 

During  the  last  several  years,  the  interest  in  fast  transforms  has  been  growing,  stimulated 
by  the  combination  of  recent  progress  in  fast  algorithms  of  various  kinds  with  the  impor¬ 
tance  of  the  Fast  Fourier  Transform  in  computational  mathematics,  electrical  engineering, 
etc.,  and  by  the  success  of  various  types  of  multilevel  computational  techniques.  In  partic¬ 
ular,  several  prior  attempts  have  been  made  to  construct  numerically  stable  fast  spherical 
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harmonic  transforms.  It  is  not  the  purpose  of  this  paper  to  provide  an  exhaustive  survey  of 
literature  on  the  subject;  we  refer  the  reader  to  the  publications  [10],  [9],  [15],  [13],  [12]  for 
details. 

The  structure  of  this  paper  is  as  follows.  In  Section  2,  we  summarize  a  number  of  facts 
(from  both  mathematical  and  numerical  analysis)  to  be  used  in  the  rest  of  the  paper;  all 
of  the  content  of  Section  2  is  either  well-known  or  follows  easily  from  well-known  facts.  In 
Section  3,  we  build  the  analytical  apparatus  to  be  used  in  the  construction  of  the  algorithms 
of  this  paper.  Section  4  contains  an  informal  description  of  the  algorithm,  and  in  Section  5 
the  procedure  is  described  in  detail.  The  performance  of  the  scheme  is  illustrated  with 
numerical  examples  in  Section  6,  and  Section  7  contains  a  discussion  of  possible  applications 
of  the  approach  of  this  paper  in  other  environments. 


2  Mathematical  and  Numerical  Preliminaries 

In  this  section,  we  summarize  several  facts  from  mathematical  and  numerical  analysis.  Please 
note  that  in  this  section  and  throughout  this  paper,  the  variable  x  always  takes  arbitrary 
values  in  [—1,  1],  9  takes  values  in  [0,  7r],  and  p  takes  values  in  [0,  2n\.  We  will  always  use 
the  term  “eigenvector”  to  mean  “normalized  eigenvector.” 

2.1  Spherical  Harmonics  and  Associated  Legendre  Functions 

In  this  subsection,  we  summarize  a  number  of  properties  of  spherical  harmonics  and  associ¬ 
ated  Legendre  functions;  all  of  these  can  be  found,  for  example,  in  [1], 

The  coefficients  in  the  spherical  harmonic  expansion  (6)  of  a  function  /  in  £2(£2)  are 
given  by  the  formula 


a 


m 

l 


bmW») 


e~imipf{9,  p)  sin  9  dp  d9. 


(13) 


For  the  forward  spherical  harmonic  transform  of  order  N,  we  have  to  compute  the  coef¬ 
ficients  (13)  from  the  values  f(9k,  ipj),  where  90 ,  9i,  ...,  9^- 1,  9N  are  defined  in  (7),  and 
Po,  ifi,  . . . ,  (f2N-i,  ^2N  are  defined  in  (8).  For  a  cost  of  0(N2  log  N),  we  use  the  Fast  Fourier 
Transform  to  obtain  the  2{2N  +  1  )(N  +  1)  values  gm{9k )  and  hm{9k )  {m  =  —N,  —IV  +  1, 
. . . ,  IV  —  1,  IV;  A;  =  0,  1,  . . . ,  IV  —  1,  IV)  of  the  functions  gm  and  hm  defined  on  [0,  ir\  by  the 


formulae 


2  TV 

9m{9)  =  ^2  cos(m  ifj)  f(9,  fj), 


3=0 


(14) 


2  N 

hm(9)  =  Y,  sin(ra  Pj)  f(9,  pj ).  (15) 

i=o 

We  then  evaluate  the  coefficients  (13)  via  the  formula 


Oil 


=  I  p\ml  (cos  9)  gm(9)  sin  9  d 9  -  i  I  P["1'  (cos  9)  hm(9)  sin  9  d9. 


H 


(16) 
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To  evaluate  the  integrals  in  (16),  we  have  to  convert  the  values  fm(cos9k)  into  the 
coefficients  p™,  p™+1,  . . . ,  p^-ij  Pn  in  expansions  of  functions  fm  on  [—1,  1]  of  the  form 

N 

rw  =  E?rprw,  (i7) 

l=m 

pT  =  j\pT(x)r(x)dx,  (18) 

where  m  is  any  integer  with  0  <  m  <  N. 

The  principal  purpose  of  this  paper  is  the  construction  of  a  “fast”  scheme  for  computing 
the  coefficients  (18)  from  the  values  /m(cos  #*,),  and  for  computing  the  inverse  of  this  trans¬ 
formation,  that  is,  for  computing  the  values  fm(cos9k)  of  the  function  fm  defined  in  (17) 
from  the  coefficients  (18). 

For  any  nonnegative  integers  l  and  m  with  m  <  l,  the  associated  Legendre  function  PJn 
on  [—1,  1]  is  defined  by  the  formula 

JTW  =  (-1  r  (19) 

where  Pi  is  the  Legendre  polynomial  of  degree  /.  Obviously,  P™  is  a  polynomial  when  m  is 
even  and  a  polynomial  multiplied  by  a/1  —  x 2  when  m  is  odd. 

For  any  nonnegative  integer  m,  we  dehne  the  differential  operator  Lm  by  the  formula 

Ufl(l)  =  ((1 " x2)  Tx!(x) ) +  rh f(x)  (20) 

for  any  function  /  on  [—1,  1]  with  a  continuous  second  derivative.  For  any  nonnegative 
integers  l  and  m  with  m  <  l ,  the  function  P™  satisfies  the  differential  equation 

Lm  (PT)  (x)  =  1(1  +  1)  P‘;\x),  (21) 

where  the  differential  operator  Lm  is  dehned  in  (20). 

For  any  integers  l  and  m  with  0  <  m  <  l  and  l  >  0, 

(2 1  +  1)  x  Pzm (p)  =  (/  +  m)  P^x)  +  (l  -  m  +  1)  P^(x).  (22) 

For  any  nonnegative  integers  l  and  m  with  1  <  m  <  l, 

/>w  <23> 

Lemma  2.1  Suppose  that  l  and  m  are  even  integers  such  that  2  <  m  <  l.  Then,  there  exist 
1/2  real  numbers  t/o,  £i,  . . . ,  C1/2- 2?  6/2-1  such  that 

1/2-1 

P?(x)  =  E  (A+2(x).  (24) 

fc= 0 
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Lemma  2.2  Suppose  that  l  and  m  are  integers  such  that  m  is  even,  l  is  odd,  and  2  <  m  <  l. 
Then,  there  exist  (/  —  l)/2  real  numbers  fo,  £i,  •  •  • ,  £(i- 1)/2-2,  £(i- i)/2-i  such  that 

__  (i-l)/2-l 

iTW=  E  (25) 

fc=0 

Lemma  2.3  Suppose  that  l  and  m  are  integers  such  that  m  is  odd,  l  is  even,  and  1  <  m  <  l. 
Then,  there  exist  1/2  real  numbers  £o,  £i,  •  •  • ,  £1/2-2,  £ 1/2-1  such  that 

1/2-1 

PT(x)  =  E  (A+2(x).  (26) 

k= 0 

Lemma  2.4  Suppose  that  l  and  m  are  odd  integers  such  that  1  <  m  <  l.  Then,  there  exist 
(/  +  l)/2  real  numbers  £0,  £ 1 ,  •  •  •,  £(i+ i)/2-2,  €(i+i)/2— 1  that 

(1+ 1)/2— 1 

E  (27) 

k=0 


2.2  Chebyshev  Polynomials 

In  this  subsection,  we  cite  the  existence  of  a  fast  algorithm  for  computing  with  Chebyshev 
polynomials. 

For  any  nonnegative  integer  k ,  we  define  Tk  to  be  the  Chebyshev  polynomial  of  degree  k 
of  the  first  kind,  defined  by  the  formula 

7fc(cos0)  =  cos  (kd)  (28) 


for  any  real  9,  and  Uk  to  be  the  Chebyshev  polynomial  of  degree  k  of  the  second  kind,  defined 
by  the  formula 


14  (cos  0) 


sin((A;  +  1)9) 
sin  9 


(29) 


for  any  real  9. 

The  following  observation  cites  the  relationship  between  the  Fast  Fourier  Transform  and 
expansions  in  series  of  Chebyshev  polynomials. 


Observation  2.5  Suppose  that  N  >  0  is  an  integer,  Co,  c\,  . . . ,  cn-i,  cn,  and  uq,  u\,  . . . , 
un-i,  un  are  real  numbers,  and  f  and  g  are  the  functions  on  [—1,  1]  defined  by  the  formulae 

N 

/(EE^W,  (30) 

k= 0 
N 

g{x)  =  E  uk  Vl  -  X2  Uk(x).  (31) 

k= 0 

Then,  there  exists  an  algorithm  which  uses  0(N  log  IV)  operations  to  convert  the  coeffi¬ 
cients  Co,  c\,  . . . ,  cn-i,  cn  into  the  values  ffx 0),  f(x  1),  . . . ,  f(xN-i),  f{xN),  and  to  convert 
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the  coefficients  u0,  U\,  . . . ,  u^-i,  ujy  into  the  values  g(x 0),  g(xf),  . . . ,  g(x n-i),  g(%N),  where 
the  sampling  locations  xq,  x\,  . . . ,  xn-i,  xn  are  defined  by  the  formula 


Xk  =  cos 


fn(k  +  1/2)  \ 

V  N  +  l  )■ 


(32) 


Moreover,  there  exists  an  algorithm  which  uses  0(N  log  N)  operations  to  convert  the  values 
f(x o),  f(x i),  . ..,  f(x n-i),  f {pc n)  into  the  coefficients  c0,  C\,  c^-i,  c n,  and  to  convert 

the  values  g(x  o),  g(x  i),  . .  .7  g(xN-i),  into  the  coefficients  uq,  u\,  . . . ,  un-i,  un  (see, 

for  example,  [If]). 


2.3  Associated  Legendre  Functions  of  Low  Orders 

In  this  subsection,  we  summarize  certain  simple  relationships  between  Chebyshev  polynomi¬ 
als  and  associated  Legendre  functions  of  orders  1  and  2.  These  relationships  are  a  straight¬ 
forward  consequence  of  formulae  7.112.1,  8.339.1,  8.339.2,  8.700.1,  8.752.1,  8.826.1,  8.828.1, 
8.832.2,  and  8.911.4  of  [7], 

We  define  the  function  A  on  [0,  oo)  by  the  formula 


AW 


r  (z  +  i) 
m+i)  ■ 


(33) 


where  Y  is  the  Euler  gamma  function. 

For  any  integer  n  >  1  and  l,  k  —  0,  1,  . . , ,  n  —  2,  n  —  1,  we  define  the  entry  A]1) ,+  of  the 
n  x  n  matrix  An'l'+  by  the  formulae 


An<  1.+ 
Al,k 


4/ +  3 


2  (2k  +  2/  +  3)  (2k  —  21  —  1) 


A(k-l)  A 


2k  +  21  +  1 


(34) 


/4(/  +  l)(2/  +  l) 
41  +  3 


when  k  >  l,  and 


A\f+  =  0 


(35) 


otherwise  (when  k  <  l). 

For  any  integer  n  >  1  and  l,  k  —  0,  1,  . . . ,  n  —  2,  n  —  1,  we  define  the  entry  A[’j]~  of  the 
n  x  n  matrix  An,l~  by  the  formulae 


An’ 1 - 
Al,k 


41  +  5 


2(2k  +  2l  +  5)(2k-2l-l) 


A  (k-l)  A 


'2k  +  21  +  3' 


(36) 


1 4(1  +  1)(2/  +  3) 
4/ +  5 


when  k  >  l,  and 
otherwise  (when  k  <  /). 


=  o 


(37) 
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For  any  integer  n  >  1  and  k,l  =  0,  1,  . . . ,  n  —  2,  n  —  1,  we  define  the  entry  B^f,+  of  the 
n  x  n  matrix  Bn,1,+  by  the  formulae 

<'+  =  2  A  (/  -  *)  A  (J  +  t  +  1)  /4(i+^+1)  (38) 

when  k  <  l,  and 

B5‘'+  =  0  (39) 

otherwise  (when  k  >  l). 

For  any  integer  n  >  1  and  k,l  =  0,  1,  . . . ,  n  —  2,  n  —  1,  we  define  the  entry  B%f’~  of  the 
n  x  n  matrix  Bn,1~  by  the  formulae 


r?n,  1,- 

Bk,i 


k  _i_  i 

4  -  A(l  —  k)  A(l  +  k  +  2) 
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I  4/  +  5 
4(/  +  1)  (21  +  3) 


(40) 


when  k  <  l,  and 

s?;,1’-  =  o  (4i) 

otherwise  (when  k  >  l). 

For  any  integer  n  >  1,  l  —  0,  1,  . . . ,  n  —  2,  n  —  1,  and  k  —  0,  1,  . . . ,  n  —  1,  n,  we  define 
the  entry  A^'£'+  of  the  n  x  (n  +  1)  matrix  ytn,2,+  by  the  formulae 


in,2,+ _  / 1  ^  2k  (6(/  +  1)(2/  +  3)  —  2(2k  —  l)(2k  +  1))  , 

^  +  (2/c  +  2/  +  3)(2/c  —  2/  —  3)  (/  j 


(42) 


A 


'2AH-2Z  +  1' 


4/ +  5 


8(/  +  !)(/  +  2)  (2 1  +  1)  (2Z  +  3) 


when  k  >  l  +  1,  and 


4*1,2,+ 

Al,k 


=  4 


4/ +  5 


8  (/  +  !)(/  +  2)(2 1  +  1)  (2Z  +  3) 


(43) 


otherwise  (when  k  <  l  +  1). 

For  any  integer  n  >  1,  /  =  0,  1,  . . . ,  n  —  2,  n  —  1,  and  k  —  0,  1,  . . . ,  n  —  1,  n,  we  define 
the  entry  of  the  n  x  (n  +  1)  matrix  An,2~  by  the  formulae 


A 


n,  2,— 
l,k 


+ 


(2k  +  1)  (6(Z  +  2) (2/  +  3)  —  8fc(fc  +  1)) 

(2k  +  2/  +  5)(2/c  —  2/  —  3)  1  J 

/  2k  +  2/  +  3\  \  I  4/  +  7 

v  2  J  J  \  8(/  +  l)(/  +  2)(2/  +  3)(2/  +  5) 


when  /c  >  l  +  1,  and 


A 


n,  2,— 
/,/e 


4/ +  7 


8(/  +  !)(/  +  2)  (2/  +  3)  (2/  +  5) 


(44) 


(45) 


otherwise  (when  &</  +  !). 


For  any  integer  n  >  1,  k  —  0,  1,  . . . ,  n  —  1,  n,  and  l  —  0,  1,  . . . ,  n  —  2,  n  —  1,  we  define 
the  entry  B^f,+  of  the  (n  +  1)  x  n  matrix  Bn’2,+  by  the  formulae 


DnA- 

Bk,i 


2(/  +  l)(2/  +  3)  -8A;2 


7 r 


A  (/  —  k  +  1)  A  (/  +  k  +  1) 


(46) 


when  k  —  0, 


]Dn,  2,- 


I  4/  +  5 

8(/  +  !)(/  +  2)(2 1  +  1)  (2Z  +  3) 


22(t  +  l)(2l  +  3)-8g 

7T 


(47) 


4(  +  5 


8(/  +  !)(/  +  2)  (2/  +  1)  (2/  +  3) 


when  0  <  k  <  l  +  1,  and 

Bnkf+  =  0  (48) 

otherwise  (when  k  >  l  +  1). 

For  any  integer  n  >  1,  k  —  0,  1,  . . . ,  n  —  1,  n,  and  /  =  0,  1,  . . . ,  n  —  2,  n  —  1,  we  define 
the  entry  B^f’~  of  the  (n  +  1)  x  n  matrix  Bn,2~  by  the  formulae 


B 


44,2,— 

k'l 


4(/  +  2)(2/  +  3)  -4(2A;  +  1)S 


7 r 


A  (/  —  fc  +  1)  A  (/  +  k  +  2) 


(49) 


4/ +  7 


8(/  +  l)(/  +  2)(2/  +  3)(2/  +  5) 


when  k  <  l  +  1,  and 

'  =  o  (50) 

otherwise  (when  /c  >  /  +  1). 

The  following  four  lemmas  are  proven  via  mechanical,  but  rather  tedious  manipulations 
of  formulae  7.112.1,  8.339.1,  8.339.2,  8.700.1,  8.752.1,  8.826.1,  8.828.1,  8.832.2,  and  8.911.4 
°f 

The  following  lemma  provides  explicit  expressions  for  the  matrix  Bn,1,+  converting  co¬ 
efficients  in  linear  combinations  of  associated  Legendre  functions  of  order  1  of  odd  degrees 
into  coefficients  in  linear  combinations  of  even  Chebyshev  polynomials  of  the  second  kind, 
scaled  by  \/l  —  x 2,  and  for  the  matrix  An,l,+  converting  the  latter  into  the  former. 


Lemma  2.6  Suppose  that  n  >  1  is  an  integer,  q  =  ( q0 ,  qi,  . . . ,  qn- 2,  qn- i)T  is  a  real  vector, 
and  f  is  the  function  on  [—1,  1]  defined  by  the  formula 

n—  1 

f(x)  =  J2(lipl2i+i(x)-  (51) 

1=0 

Then, 

n—  1 

f(x)  =  Y1  ukVl-  X 2  U2k(x ),  (52) 

k= 0 
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(53) 


where  u  =  ( u0 ,  ult  . . . ,  un_ 2,  wn-i)T  is  the  real  vector  defined  by  the  formula 

u  =  Bn,1,+  q, 

and  Bn,1,+  is  defined  in  (33),  (38),  and  (39).  Furthermore, 

q  =  AnX+u ,  (54) 

where  AU)l'+  is  defined  in  (33),  (34),  mid  (35). 

The  following  lemma  provides  explicit  expressions  for  the  matrix  B n,1~  converting  co¬ 
efficients  in  linear  combinations  of  associated  Legendre  functions  of  order  1  of  even  degrees 
into  coefficients  in  linear  combinations  of  odd  Chebyshev  polynomials  of  the  second  kind, 
scaled  by  a/1  —  x 2,  and  for  the  matrix  An’l,~  converting  the  latter  into  the  former. 

Lemma  2.7  Suppose  that  n  >  1  is  an  integer,  q  =  ( q0 ,  qi,  . . . ,  qn-2,  Qn- i)T  is  a  real  vector, 
and  f  is  the  function  on  [—1,  1]  defined  by  the  formula 

n— 1 

f{x)  =  YJ(liP\i+2^)-  (55) 

1=0 

Then, 

n— 1 

f(x)  =  Y,  ukVl-  X 2  U2k+i(x),  (56) 

k= 0 

where  u  =  ( u0 ,  ult  . . . ,  un_ 2,  un-i)T  is  the  real  vector  defined  by  the  formula 

u  =  Bn,1~  q,  (57) 

and  Bn,1~  is  defined  in  (33),  (40),  and  (41)-  Furthermore, 

q  =  An'l~  u,  (58) 

where  An,1~  is  defined  in  (33),  (36),  and  (37). 

The  following  lemma  provides  explicit  expressions  for  the  matrix  Bn,2,+  converting  co¬ 
efficients  in  linear  combinations  of  associated  Legendre  functions  of  order  2  of  even  degrees 
into  coefficients  in  linear  combinations  of  even  Chebyshev  polynomials  of  the  first  kind,  and 
for  the  matrix  An,2,+  converting  the  latter  into  the  former. 

Lemma  2.8  Suppose  that  n  >  1  is  an  integer,  p  =  ( p0 ,  pi,  . . . ,  pn_ 2,  pn- i)T  is  a  real  vector, 
and  f  is  the  function  on  [—1,  1]  defined  by  the  formula 

n—  1 

f(%)  =  YPiph+ 2  0*0-  (59) 

1=0 

Then, 

n 

f(X)  =  YCkT2k{x),  (60) 

k= 0 
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where  c  =  (c0,  Ci,  . . . ,  cn_i,  cn)T  is  t/ie  rea/  vector  defined  by  the  formula 

c  =  Bn'2’+p ,  (61) 

and  Bn,2,+  is  defined  in  (33),  (46),  (4V>  am d  (4$)-  Furthermore, 

p  =  An’2’+  c,  (62) 

where  An,2’+  is  defined  in  (33),  (42),  and  (43). 

The  following  lemma  provides  explicit  expressions  for  the  matrix  Bn'2~  converting  co¬ 
efficients  in  linear  combinations  of  associated  Legendre  functions  of  order  2  of  odd  degrees 
into  coefficients  in  linear  combinations  of  odd  Chebyshev  polynomials  of  the  first  kind,  and 
for  the  matrix  An,2~  converting  the  latter  into  the  former. 

Lemma  2.9  Suppose  that  n  >  1  is  an  integer,  p  =  (po,  pi,  . . . ,  pn- 2,  Pn- i)T  is  a  real  vector, 
and  f  is  the  function  on  [—1,  1]  defined  by  the  formula 

n—  1 

f(x)  =  YjPiP221+Ax)-  (63) 

1=0 

Then, 

n 

f(x)  =  Y,ck  T2k+i(x ),  (64) 

k= 0 

where  c  =  (c0,  Ci,  . . . ,  cn_  1,  cn)T  is  the  real  vector  defined  by  the  formula 

c  =  Bn,2~  p,  (65) 

and  Bn’2~  is  defined  in  (33),  (49),  and  (50).  Furthermore, 

p  =  An’2’~  c,  (66) 

where  An,2~  is  defined  in  (33),  (44):  and  (45)- 

Observation  2.10  Suppose  that  n  >  1  is  an  integer.  Then,  there  exists  an  algorithm  which 
uses  0(n  log (n/e))  operations  to  apply  to  an  arbitrary  vector  any  of  the  matrices  An'1,+ , 
BnF+ ,  A11'1- ,  BnF~ ;  An’2,+,  Bn’2>+,  An’2’~,  and  Bn’2~  defined  in  (33)-(50),  where  £  is  the 
precision  of  computations  (see  [3]). 


2.4  Semiseparable  Matrices 

For  any  integer  n  >  0,  a  semiseparable  real  n  x  n  matrix  S'  is  a  matrix  whose  entry  Sj^  is 
given  by  the  formulae 

Sj,k  =  dj  bk  (67) 


when  j  <  k,  and 


Fj.k  bj 


(68) 
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when  j  >  k,  where  a  =  ( a0 ,  a\,  . . . ,  an_2,  an-i)T  and  b  =  ( b0 ,  bi,  . . . ,  bn_2,  &n-i)T  are  real 
vectors. 

Matrices  of  the  form 

G  =  D  +  S,  (69) 

where  D  is  a  diagonal  real  matrix  and  S'  is  a  semiseparable  real  matrix,  will  be  encountered 
repeatedly  throughout  this  paper.  The  matrix  U  of  eigenvectors  of  the  matrix  G  in  (69)  will 
be  particularly  important;  U  is  orthogonal  and  diagonalizes  G,  so  that 

UT  GU  =  A,  (70) 


where  A  is  a  diagonal  real  matrix. 

The  principal  numerical  tool  of  this  paper  is  the  following  observation,  made  in  [8]  and  [4] . 

Observation  2.11  The  matrices  U  and  UT  in  (70)  can  be  applied  to  an  arbitrary  vector  of 
length  N  for  a  cost  of  0(N (log  N)  log(l/e))  operations,  where  e  is  the  precision  of  compu¬ 
tations. 

Remark  2.12  Strictly  speaking,  only  the  numerical  apparatus  behind  Observation  2.11 
is  constructed  in  [4],  However,  the  observation  itself  is  stated  explicitly  in  a  very  similar 
environment  in  [8].  In  our  implementation,  we  used  a  minor  modification  of  the  apparatus 
in  [4],  to  be  reported  at  a  later  date. 

3  Analytical  Apparatus 

In  this  section,  we  construct  the  principal  analytical  tools  used  in  this  paper. 

In  Subsection  3.1,  we  observe  that  when  the  function  P™  is  represented  as  a  linear 
combination  of  functions  P-  or  P)  (depending  on  whether  m  is  even  or  odd),  the  Sturm- 
Liouvillc  problem  (21)  becomes  an  eigenvector  problem  for  the  matrix  G  in  (69).  Thus, 
according  to  Observation  2.11,  there  exists  an  algorithm  that  uses  O  (iV(log  N)  log(l/e)) 
operations  to  apply  the  matrices  U  and  UT  in  (70)  to  arbitrary  vectors  of  length  N ,  where 
£  is  the  precision  of  computations. 

In  Subsection  3.2,  we  observe  that  the  problem  of  evaluating  expansions  of  the  form  (12) 
can  be  reduced  to  the  problem  of  evaluating  expansions  of  the  forms  (10)  and  (11),  via  the 
matrices  U  and  UT  in  (70).  These  matrices  can  be  applied  to  arbitrary  vectors  efficiently, 
due  to  Observation  2.11. 

3.1  Associated  Legendre  Differential  Equations  in  Terms  of  Asso¬ 
ciated  Legendre  Functions  of  Low  Orders 

For  any  even  integers  n  and  m  with  2  <  m  <  n,  and  for  j,k  =  0,  1,  . . . ,  n / 2  —  2,  n/2  —  1, 
we  define  the  entry  Gn-'ffl  of  the  n/2  x  n/2  matrix  Gn,m  by  the  formula 

G”f  =  £  Ff+2(x)  Lm  (P22t+2)  (x)  dx,  (71) 
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where  the  differential  operator  Lm  is  defined  in  (20). 

For  any  odd  integer  n  and  even  integer  m  with  2  <  m  <  n,  and  for  j,k  =  0,  1,  . . . , 
(n  —  l)/2  —  2,  (n  —  l)/2  —  1,  we  dehne  the  entry  G£™  of  the  (n  —  l)/2  x  (n  —  l)/2  matrix 

Qn,m  ^ie  formuia 

G”f  =  £  F%+  3(x)  Lm  ( P22t+3 )  (x)  dx,  (72) 

where  the  differential  operator  Lm  is  dehned  in  (20). 

For  any  even  integer  n  and  odd  integer  m  with  1  <  m  <  n,  and  for  j,k  =  0,  1,  . . . , 
nj 2  —  2,  n/2  —  1,  we  dehne  the  entry  G”’™  of  the  n/2  x  nj 2  matrix  G",m  by  the  formula 

G”f  =  £  P\s+2(x)  Lm  (P‘t+2)  (x)  dx ,  (73) 

where  the  differential  operator  Lm  is  dehned  in  (20). 

For  any  odd  integers  n  and  m  with  1  <  m  <  n,  and  for  j,  k  —  0,  1,  . . . ,  (n  +  l)/2  — 
2,  (n  +  l)/2  —  1,  we  dehne  the  entry  G”’£l  of  the  (n  +  l)/2  x  (n  +  l)/2  matrix  Gn,m  by  the 
formula 

G”f  =  £  p\j+1(x)  Lm  (PGi)  (x)  dx,  (74) 

where  the  differential  operator  Lm  is  dehned  in  (20). 

The  following  lemma  states  that  the  coefficients  in  the  expansion  of  the  function  P£  in 
terms  of  either  the  functions  P  ■  or  the  functions  P  •  (depending  on  whether  m  is  even  or 
odd)  are  the  entries  in  an  eigenvector  of  the  matrix  Gn,m. 

Lemma  3.1  Suppose  that  m  and  n  are  integers  such  that  1  <  m  <  n. 

Then,  when  m  and  n  are  both  even,  1(1  +  1)  is  an  eigenvalue  of  the  matrix  Gn,m  defined 
in  (71 ),  for  any  even  integer  l  with  m  <  l  <  n,  and  the  coordinates  £0,  £i,  . . . ,  fn/2-2 ,  fn/ 2-1 
of  the  corresponding  eigenvector  are  the  coefficients  in  the  expansion 

n/ 2-1 

PT(X)=  Y.  (A+2(x),  (75) 

k= 0 

When  m  is  even  and  n  is  odd,  1(1  +  1 )  is  an  eigenvalue  of  the  matrix  Gn,m  defined  in  (72), 
for  any  odd  integer  l  with  m  <  l  <  n,  and  the  coordinates  £0,  £1,  . . .  7  £>-i)/2-i 

of  the  corresponding  eigenvector  are  the  coefficients  in  the  expansion 

(n- 1)/2— 1 

P?(x)=  E  (W2t+3(x).  (76) 

k= 0 

When  m  is  odd  and  n  is  even,  1(1  +  1 )  is  an  eigenvalue  of  the  matrix  Gn,m  defined  in  (73), 
for  any  even  integer  l  with  m  <  l  <  n,  and  the  coordinates  £).,  £1,  •  •  • ..  in/2-2,  £7/2-1  of  the 
corresponding  eigenvector  are  the  coefficients  in  the  expansion 

n/ 2-1 

P)(x)=Y.i*P+ ,0).  (77) 

k= 0 
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When  m  and  n  are  both  odd,  l{l  +  1)  is  an  eigenvalue  of  the  matrix  Gn,m  defined  in  (74), 
for  any  odd  integer  l  with  m  <1  <  n,  and  the  coordinates  6n  6l>  . . . ;  £(n+i)/2-2,  £(n+ i)/2-i 
of  the  corresponding  eigenvector  are  the  coefficients  in  the  expansion 

(n+l)/2— 1 

P?0)=  E  (78) 

k= 0 

Proof.  We  outline  the  proof  in  the  case  that  m  and  n  are  both  even;  the  proofs  in  the  other 
three  cases  are  similar. 

Substituting  (24)  into  (21)  and  using  (2)  and  (5),  we  obtain  from  (71)  that  the  numbers 
6),  fi,  . . . ,  6/2-2,  6/2-1  from  (24),  along  with  the  numbers  6/2  =  0,  6/2+1  =  0,  . . . ,  fn/2-2  = 
0,  61/2—1  =  0  when  l  <  n,  are  the  coordinates  of  the  eigenvector  of  Gn,m  with  corresponding 
eigenvalue  /(/  +  1),  giving  the  expansion  (75).  □ 

The  following  lemma  states  that  Gn'm  is  the  sum  of  a  diagonal  matrix  and  a  semiseparable 
matrix  and  provides  expressions  for  the  entries  of  Gn,m. 

Lemma  3.2  Suppose  that  m  and  n  are  integers  such  that  1  <  m  <  n. 

Then,  the  matrix  Gn,m  defined  in  (71)-(74)  has  the  form 

Gn'm  =  D  +  S,  (79) 

where  D  is  a  diagonal  matrix  and  S  is  a  semiseparable  matrix. 

When  m  and  n  are  both  even,  D  is  the  diagonal  n/ 2  x  n j 2  matrix  with  the  diagonal 
entries  D0fi,  Dhl,  ...,  Dn/2_2^/2_2,  Dn/2-i,n/2- 1  defined  by  the  formula 

Dk:k  =  (2k  +  2)(2k  +  3),  (80) 

and  S  is  the  semiseparable  n/2  x  n/2  matrix  with  the  entry  Sj:k  defined  by  the  formulae 

Sj,k  =  dj  bk  (81) 

when  j  <  k,  and 

Sj)k  dk  bj 

otherwise  (when  j  >  k),  where  the  numbers  a0,  0 1 ,  ...,  an/2_ 2,  an/ 2-1  and  b0,  b\, 
bn/2-2,  bn/ 2-1  are  defined  by  the  formulae 


ak  — 


(2k  +  1)  (27c  +  2)(2  k  +  3)(2fc  +  4)(4A:  +  5) 


8  •  15 


(83) 


bk  =  (' m 2  -  4) 


15(4*;  + 5) 


(84) 


2(2k  +  1)(2 k  +  2)(2 k  +  3) (2 k  +  4) ' 

When  m  is  even  and  n  is  odd,  D  is  the  diagonal  (n  —  l)/2  x  (n  —  l)/2  matrix  with 
the  diagonal  entries  Dq  q,  D\  \,  . . . ,  D^n_ i)/2— 2,(n— 1)/2— 2,  D(n— 1)/2— i,(n— 1)/2— 1  defined  by  the 
formula 

Dk,k  =  (2k  +  3)(2k  +  A),  (85) 
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and  S  is  the  semiseparable  (■ n  —  l)/2  x  (n  —  l)/2  matrix  with  the  entry  Shk  defined  by  the 
formulae 

aj  bk  (86) 

when  j  <  k, 

Sj,k  =  ak  bj  (87) 

otherwise  (when  j  >  k),  where  the  numbers  a0,  cii,  . . .  7  a(n_i)/2_2,  a(n-i)/2-i  and  b0,  bi,  . . . , 
b(n- 1)/2_2,  6(n- i)/2-i  are  defined  by  the  formulae 


dk 


(2k  +  2)  (2  k  +  3)  (27c  +  4)  (2  k  +  5)(4  k  +  7) 


14  •  15 


(88) 


bk  =  (m2  -  4) 


(89) 


7  •  15(4ife  +  7) 

AJ  8(2fc  +  2) (2k  +  3)(2fc  +  4)(2fc  +  5) ' 

When  m  is  odd  and  n  is  even,  D  is  the  diagonal  n/2  x  nj 2  matrix  with  the  diagonal 
entries  D0fi,  Dhl,  ...,  -D„/2-2,n/2-2,  Dn/2- i,„/2-i  defined  by  the  formula 

Dk>k  =  (2k  +  2)(2k  +  3),  (90) 

and  S  is  the  semiseparable  n/2  x  n/2  matrix  with  the  entry  S])k  defined  by  the  formulae 

Sj,k  =  Oj  bk  (91) 

when  j  <  k,  and 


. /,■  ak  bj 

otherwise  (when  j  >  k),  where  the  numbers  ao,  oi 
bn/2-2,  bn/ 2-1  are  defined  by  the  formulae 


?  *  *  *  ; 


an/2-2,  an/2-i  and  bo,  b\ 


ak  — 


(2k  +  2)  (27c  +  3)(4/c  +  5) 
30  : 


=  (m2  -  1) 


15(4/c  +  5) 


(92) 


?  *  *  *  ; 


(93) 


(94) 


\|  2(2fc  +  2)(2fc  +  3) ' 

When  m  and  n  are  both  odd,  D  is  the  diagonal  (n  +  l)/2  x  (n  +  l)/2  matrix  with  the 
diagonal  entries  Dog,  D\ik,  ■  ■  ■ ,  D{n+i)/2— 2,(n+i)/2— 2,  D(n+ 1)/2— i,(n+i)/2— 1  defined  by  the  for¬ 
mula 

Dk>k  =  (2k  +  l)(2k  +  2),  (95) 

and  S  is  the  semiseparable  (n  +  l)/2  x  (n  +  l)/2  matrix  with  the  entry  Sj ^  defined  by  the 
formulae 

Hj . /,■  dj  bk  (96) 

when  j  <  k,  and 

Sj,k  =  ak  bj  (97) 
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otherwise  (when  j  >  k),  where  the  numbers  a0,  a\,  . . . ,  a(n+i)/2-2>  a{n+ i)/2-i  and,  b0,  b\,  . . . , 
b(n+ 1)/2-2,  6(n+i)/2-i  ore  defined  by  the  formulae 


(2k  +  l)(2/c  +  2)(4/c  +  3) 
6  : 


bk  =  (m2  -  1) 


3(4k  +  3) 

\|  2(2/c  +  l)(2k  +  2) ' 


(98) 

(99) 


Proof.  We  outline  the  proof  in  the  case  that  m  and  n  are  both  even;  the  proofs  in  the  other 
three  cases  are  similar. 

We  define  the  entry  Dj  k  of  the  n / 2  x  n/2  matrix  D  by  the  formula 


Djtk  — 


P 


i-i 


2j+2 


X 


(p 


2 

2k+2 


[X\ 


dx , 


(100) 


where  the  differential  operator  L2  is  defined  in  (20),  and  we  define  the  entry  Sj  k  of  the 
n/2  x  n/2  matrix  S  by  the  formula 


Sj,k  — 


m2  —  4 
1  —  x2 


P22t+2(x)dx. 


(101) 


We  now  show  that  Gn,m  =  D  +  S,  D  is  diagonal,  and  S  is  semiseparable. 

Combining  (71),  (100),  (101),  and  (20),  we  obtain  the  decomposition  (79). 

Substituting  (21)  into  (100),  and  using  (2)  and  (5),  we  observe  that  the  matrix  D  is 
diagonal,  with  the  diagonal  entries  given  by  (80). 

In  order  to  obtain  the  formulae  (81)-(84),  we  define  the  entry  M])k  of  the  infinite¬ 
dimensional  matrix  M,  for  j,k  =  0,  1,2,  . . . ,  by  the  formula 

/I  _ 9  —  4  _ 9 

P-2j+2(x)  _  2  P2k+ 2O)  dx,  (102) 

-  1  -L  Jb 


and  observe  that  the  entry  (M  1)jyk  of  the  inverse  M  1  of  the  matrix  M  is  given  by  the 
formula 

(M-\t  =  //  P22j+2(x )  p\t+2(x)  dx,  (103) 

since  M  represents  the  operator  acting  on  functions  on  [—1,  1]  by  multiplication  by  the  factor 


m2  —  4 
1  —  x2  ’ 


(104) 


whereas  M  1 
inverse  factor 


represents  the  operator  acting  on  functions  on  [—1,  1]  by  multiplication  by  the 

1  —  x2 
m2  —  4 


Using  (2),  (5),  (22),  and  (103),  we  observe  that  M  1  is  tridiagonal.  So,  M  is  the  inverse 
of  a  tridiagonal  matrix,  and,  as  such,  M  is  semiseparable  (see,  for  example,  [6]).  But,  for 
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j,k  =  0,  1,  . . . ,  n / 2  —  2,  n/2  —  1,  (101)  and  (102)  show  that  Sj)k  =  M]k,  so  that  S  is  also 
semiseparable. 

Integrating  by  parts  a  few  times,  while  using  (2),  (5),  and  (19),  together  with  (101), 
yields  explicit  expressions  for  the  leftmost  entries  So,o,  *So,i,  . . . ,  50>n/ 2-2,  <S'o,n/2-i  °f  S.  Com¬ 
bining  (2),  (5),  (23),  and  (101)  yields  explicit  expressions  for  the  diagonal  entries  So,o,  <Si,i, 
. . . ,  Sn/2—2,n/2—2 ;  <S'n/2-i,n/2-i  °f  S.  Combining  all  of  these  explicit  expressions  with  the  fact 
that  S  is  semiseparable,  we  obtain  the  formulae  (81)-(84).  □ 


3.2  Associated  Legendre  Expansions  of  Arbitrary  Orders  and  As¬ 
sociated  Legendre  Expansions  of  Low  Orders 

Lemmas  3.3,  3.4,  3.5,  and  3.6  of  this  subsection  follow  immediately  from  Lemma  3.1. 

Lemma  3.3  states  that  the  matrix  of  eigenvectors  of  the  matrix  Gn,m  in  (79)  converts 
the  coefficients  in  the  expansion  of  a  function  /  in  terms  of  the  functions  P™ ,  P™+2,  . . . , 

P™_ 2,  P™  into  the  coefficients  in  the  expansion  of  the  function  /  in  terms  of  the  functions 
— 2  — 2  — 2  — 2 

P2,  P4,  . . . ,  Pn_ 2,  Pn.  Lemma  3.3  states,  moreover,  that  the  adjoint  of  the  matrix  of 
eigenvectors  converts  the  latter  coefficients  into  the  former  coefficients. 


Lemma  3.3  Suppose  that  n  and  m  are  even  integers  such  that  2  <  m  <  n,  pm  =  (p™,  p™, 
■■■>  Pn/2-21  Pn/ 2-i)T  is  a  real  column  vector  such  that  p^_m)/2+l  =  0,  p™n_m)/2+2  =  0,  . . . , 
Prn/2-2  =  Pn/ 2-1  =  0?  cind  f  is  the  function  defined  on  [—1,  1]  by  the  formula 

(n—m)/ 2 

f(x)=  Y,  P?P2j+rniX)-  (106) 

3=0 

Then, 

n/ 2-1 

f[x)  =  X]  Plph+ 2O),  (107) 

1=0 

where  p2  =  (p%,  p\,  . . . ,  p2n/2_2 ,  p2n/ 2_i)t  is  the  real  vector  defined  by  the  formula 

p2  =  Upm,  (108) 

and  U  is  an  n/2  x  n/2  matrix  of  eigenvectors  of  the  symmetric  matrix  Gn'm  in  (79),  with 

ra/2-1 

K+m(x)=  £  U,jPll+2(x)  (109) 

1=0 

for  j  =  0,  1,  . . . ,  (n  —  m)/2  —  1,  (n  —  m)/2.  Moreover, 

pm  =  UTp2.  (110) 


Lemmas  3.4,  3.5,  and  3.6  are  analogues  of  Lemma  3.3  for  different  conditions  on  m  and 


n. 
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Lemma  3.4  Suppose  that  n  and  m  are  integers  such  that  m  is  even,  n  is  odd,  2  <  m  < 
n,  pm  =  (Po*,  pT,  •••;  V™n- 5)/2>  P(^_3)/2)T  a  real  cofamn  rector  such  that  p™n_m+l)/2  = 
0,  p^_m+3)/2  =  0,  . . .  7  p^_5)/2  =  0,  pf/_3)/2  =  0,  and  f  is  the  function  defined  on  [- 1,  1]  by 
the  formula 


(n—m— 1)/2 

/(*)  =  £  p”n”+m+.M- 

j-o 


Then, 


where  p2  =  (pg,  Pi, 


(n- 3)/2 


/<X>  =  ^  P^i+3^)) 


(111) 


(112) 


«=o 


v\n- 5)/2)  P(„._3)/2)T  rea^  vector  defined  by  the  formula 

p2  =  Upm ,  (113) 


and  U  is  an  (■ n  —  l)/2  x  (n  —  l)/2  matrix  of  eigenvectors  of  the  symmetric  matrix  Gr 
in  (79),  with 

(n— 3)/2 


P 


0=  £ 

J=0 


2j+m+l 

/or  /  =  0,  1,  . . . ,  (n  —  m  —  3)/2;  (n  —  m  —  l)/2.  Moreover. 

pm  =  jjT  p2 ' 


X) 


(114) 


(115) 


Lemma  3.5  Suppose  that  n  and  m  are  integers  such  that  m  is  odd,  n  is  even,  1  <  m  < 

n,  pm  —  (p™,  p™,  . ..,  P™/2_2,  P™/2_i)t  a  real  column  vector  such  that  pZ_m+iy2  = 

o,  P(n— m+3)/2  =  °>  •••;  tC/2-2  =  0,  P™/i-i  =  0,  and  f  is  the  function  defined  on  [-1,  1] 
by  the  formula 


(n—m— 1)/2 

/(*)  =  £  p”^+m+,W. 

3-0 


(116) 


Then, 


n/ 2-1 

/(a;)  =  ^  P;1?^^).  (117) 

1=0 

where  p1  =  (pj,  pj,  . . . ,  p^/2_2,  Pn/2-i)T  zs  t/ie  rea/  vector  defined  by  the  formula 

p1  =  f/pm,  (118) 

and  U  is  an  nj 2  x  n / 2  matrix  of  eigenvectors  of  the  symmetric  matrix  Gn,m  in  (79),  with 

n/ 2-1 

(119) 


1=0 


x) 


for  j  =  0,  1,  . . . ;  (n  —  m  —  3)/2;  (n  —  rn  —  l)/2.  Moreover, 


pm  =  UTp\ 


(120) 
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Lemma  3.6  Suppose  that  n  and  m  are  odd  integers  such  that  1  <  m  <  n,  pm  =  (p™,  pfi, 
P™n- 3)/2.  Pfn-1)/2)T  is  a  real  column  vector  such  that  p™n_m)/2+1  =  0,  p™n_m)/2+2  =  0,  . . . , 
P{n- 3)/2  =  0’  Pon- 1)/2  =  0>  an(h  /  the  function  defined  on  [—1,  1]  by  the  formula 


(n—m)/2 

f(X)=  Y  PT^j+m(X)- 

3= 0 


Then, 


(n-l)/2 


/(»=  Z  Pj-Pa+iW. 


(121) 


(122) 


«=o 


where  p1  =  (pg,  p*,  . . .  ?  P(n_3\/2,  P(n-i)/2)T  the  real  vector  defined  by  the  formula 

p1  =  Upm ,  (123) 

and  U  is  an  (n  +  l)/2  x  (n  +  l)/2  matrix  of  eigenvectors  of  the  symmetric  matrix  Gn,m 
in  (79),  with 

(n— 1)/2 


Py+min  =  £ 

z=o 

for  j  =  0,  1,  . . . ,  (n  —  m)/ 2  —  1,  (n  —  m)/ 2.  Moreover, 

pm  =  nTpi 


x 


(124) 


(125) 


4  Informal  Description  of  the  Algorithm 


In  this  section,  we  outline  a  “fast”  algorithm  for  the  conversion  of  the  values  /g,  fi,  . .., 
/jv-i,  /v  °f  the  function  /  tabulated  at  the  nodes  a;0,  x±,  . ..,  xjv-i,  atw,  defined  by  the 
formula 


Xk  =  cos 


/  ?r(A:  +  1/2)  \ 

V  N  +  1  )' 


(126) 


into  the  coefficients  Po,  Pi,  ■  ■  ■ ,  PN-m- i,  PN-m  hr  the  expansion 


N—m 

}{X)=  Y  PlPT+m(X)- 


1=0 


(127) 


(The  inverse  procedure  of  converting  the  coefficients  po,  pi,  . . . ,  pN-m-i,  PN-m  into  the 
values  /o,  fi,  _ ,  /v-i,  /at  is  quite  similar,  so  we  omit  its  description.) 

The  procedure  consists  of  four  steps,  described  briefly  below.  Since  the  procedure  when 
m  is  odd  is  virtually  the  same  as  when  m  is  even,  we  describe  the  procedure  only  for  the 
case  that  m  is  even.  For  definitiveness,  we  assume  also  that  N  is  even. 


1.  We  separate  /  into  its  even  and  odd  parts  f+  and  /  ,  defined  by  the  formulae 

f(x)  +  f{~x) 


f+(x)  = 


V2 


(128) 
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(129) 


/  (a)  = 


f(x)  -  f(~x) 

V2 


All  subsequent  processing  is  performed  separately  for  f+  and  /  .  Since  the  procedures  for 
f+  and  f~  are  virtually  identical,  we  describe  only  the  procedure  for  the  even  part  f+. 


2.  Using  the  Fast  Fourier  Transform  as  in  Observation  2.5,  we  convert  the  values  /q",  /{*", 
. . . ,  /tv_i,  In  °f  the  function  /+,  defined  in  (128),  at  the  points  xo,  xi,  . . . ,  xn- i,  xn,  defined 
in  (126),  into  the  coefficients  c0,  Ci,  . . . ,  cjv/2-i)  cn/2  in  the  expansion 

N/2 

f+{x )  =  J2Ck  T'2k(x).  (130) 

k= 0 


3.  Using  Lemma  2.8  and  Observation  2.10,  we  convert  the  coefficients  c0,  Ci,  . . . ,  cjv/2-1,  c^/2 
in  the  expansion  (130)  into  the  coefficients  p%,  p\,  . . . ,  p2N/ 2_2,  p2N/2-\  in  the  expansion 

N/2—1 

f+(x)=  E  Pl?2!+2(4  (131) 

1=0 


4.  Using  Lemmas  3.2  and  3.3  and  Observation  2.11,  we  convert  the  coefficients  pf},  p\,  . . . , 
Pn/ 2-2)  Pn/ 2-1  in  the  expansion  (131)  into  the  coefficients  p £\  pf,  . . . ,  ,  P(w-m)/ 2 

in  the  expansion 

(N—m)/2 

f\x)=  Y,  pTP S+JU  (132) 

1=0 

Remark  4.1  Clearly,  Step  1  takes  O(N)  operations.  As  per  Observation  2.5,  Step  2 
takes  0(N  log  At)  operations.  As  per  Observation  2.10,  Step  3  takes  0(N  log  (At/e))  op¬ 
erations,  where  e  is  the  precision  of  computations.  As  per  Observation  2.11,  Step  4  takes 
d?(AT(log  AT)  log(l/e))  operations,  where  e  is  the  precision  of  computations  used  in  this  step. 
All  together,  Steps  1-4  take  C?(At(log  At)  log(l/e))  operations. 

Remark  4.2  To  handle  f~  when  m  is  even,  we  substitute  Lemma  2.9  for  Lemma  2.8  in 
Step  3  and  Lemma  3.4  for  Lemma  3.3  in  Step  4.  For  f~  when  m  is  even,  we  find  that 
Pn/2—1  =  0  and  cn/2  =  0,  since,  when  m  is  even,  f~  is  a  linear  combination  of  only  N/2  —  1 

functions  P3,  P5,  . . . ,  PN_ 3,  U7V_1  (or,  equivalently,  of  the  N/2  Chebyshev  polynomials  of 
the  first  kind  Tu  T3,  . . . ,  TN_ 3,  TN_/). 

To  handle  f+  when  m  is  odd,  we  substitute  Lemma  2.6  for  Lemma  2.8  in  Step  3  and 
Lemma  3.5  for  Lemma  3.3  in  Step  4. 

To  handle  f~  when  m  is  odd,  we  substitute  Lemma  2.7  for  Lemma  2.8  in  Step  3  and 
Lemma  3.6  for  Lemma  3.3  in  Step  4. 

When  computing  the  coefficients  (13)  in  the  spherical  harmonic  expansion  (6),  we  tab¬ 
ulate  the  function  /  on  [—1,  1]  at  the  same  N  +  1  sampling  locations  Xo,  Xi,  . . . ,  Xn-i,  Xn 
when  m  is  odd  as  when  m  is  even,  even  though  when  m  is  odd,  /  is  a  linear  combination  of 
only  N  functions  P3,  P2,  ■  ■  ■ ,  Pn-ii  P n  (°U  equivalently,  of  the  N  Chebyshev  polynomials 
of  the  second  kind  y/l  —  x2Uo,  \/l  —  x2  U\,  ...,  \/l  —  x2  Un-2-,  \/l  —  x2  Un-i,  which  are 
scaled  by  \/l  —  x2). 
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5  Detailed  Description  of  the  Algorithm 

In  this  section,  we  describe  in  detail  the  algorithm  described  informally  in  Section  4. 

Precomputations. 

Comment  [Compute  all  the  data  that  Steps  3  and  4  require  that  do  not  depend  on  /.] 
Compute  all  the  data  that  Step  3  requires  to  apply  fast  the  matrix  AN^2,2,+  from  Lemma  2.8 
(he.,  “compress”  the  matrix  AN/2,2,+). 

Compute  all  the  data  that  Step  4  requires  to  apply  fast  the  adjoint  of  the  matrix  U  from 
Lemmas  3.2  and  3.3  (he.,  “compress”  the  matrix  UT). 

Step  1. 

Comment  [Convert  the  values  /0,  A,  •  •  • ,  /jv-i,  /tv  into  the  values  /0+,  ./j+,  . . . ,  /^_1;  f£.] 
do  n  =  0,  . . . ,  N 

Set  fn  =  ( fn  +  fN-n)fV 2- 

enddo 


Step  2. 

Comment  [Convert  the  values  /q1",  // ,  ...,  //_l5  //  into  the  coefficients  Co,  Ci,  ..., 
Ctv/2— 1  j  C7V/2  of  the  Chebyshev  polynomials  T0,  T2,  . . . ,  Tv-2,  Tv-] 

Use  the  Discrete  Cosine  Transform  (see,  for  example,  [14])  to  convert  the  values  //,  // ,  . . . , 
Tv-i,  /tv  into  tlie  coefficients  c0,  ci,  . . . ,  cat/2-i,  cat/2- 


Step  3. 

Comment  [Convert  the  coefficients  Cq,  ci,  . . . ,  Cjv/2-1,  Cat/2  of  the  Chebyshev  polynomials 
T0,  T2,  . . . ,  Tv-2,  Tv  into  the  coefficients  pjj,  Pi,  •  •  • ,  Ptv/2-2,  Pn/2-i  °f  the  associated  Legen¬ 
dre  functions  P2,  -P4,  •  •  • ,  Pjv-2,  -^tv-] 

Use  Observation  2.10  to  apply  the  matrix  AN/2,2,+  from  Lemma  2.8  to  the  vector  c  = 
(c0,  ci,  ... ,  Cjv/2— 1,  cjv/2)T,  in  order  to  obtain  the  vector  p2  =  (p%,  pj,  . . . ,  p^/2_2,  pir/2_i)T- 

Step  4. 

Comment  [Convert  the  coefficients  Pq,  Pi,  •  •  • ,  Pat/2-2,  Pat/ 2-1  °f  the  associated  Legendre 
functions  P2,  P4,  . . . ,  P^_2,  P^  into  the  coefficients  p£\  p™,  . . . ,  P^_m)/2_i,  Pfv-m)/2  of 
the  associated  Legendre  functions  P™,  P™+2 ,  . . . ,  Pjy_2,  P^.] 

Use  Observation  2.11  to  apply  the  adjoint  of  the  matrix  U  from  Lemmas  3.2  and  3.3  to  the 
vector  p 2  =  (Pq,  Pi,  •  •  • ,  Pat/2-2,  Ptv/2-i)T  in  order  to  obtain  the  vector  pm  =  (p™,  pjn,  . . . , 

\T 

P TV/2— 2,  PtV/2— 1/  ’ 

Comment  [Only  the  first  (AT  —  m)/2  +  1  entries  of  pm  interest  us;  the  other  entries  vanish: 
P(AT— m)/2+l  =  0,  P(AT-m)/2+2  =  0,  '  '  '  ,  PtV/2-2  =  PjV/2-1  =  O'] 
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6  Numerical  Results 


The  algorithms  described  in  this  paper  have  been  implemented  in  Fortran.  Tables  1-8  below 
report  the  results  of  applying  the  algorithm  to  functions  /  defined  on  [—1,  1]  by  the  formulae 

(N—m) /2 

/M=  £  (133) 

1=0 

when  m  is  even  and  degrees  are  even, 

(N—m— 2) /2 

f(x)  =  Y  £lPll+m+l(X)  (134) 

1=0 

when  m  is  even  and  degrees  are  odd, 

(JV— m— 1)/2 

f(x)=  Y  £lPZl+m+l(X)  (135) 

1=0 

when  m  is  odd  and  degrees  are  even,  and 

(AT— m— 1)/2 

f(x)  =  Y  £lP7l+m(X)  (136) 

1=0 

when  m  is  odd  and  degrees  are  odd;  Eq,  e i,  £2,  ...  are  randomly  generated  numbers  in  the 
interval  [—1,  1].  The  decomposition  algorithm  computes  from  the  sample  values  of  /  the 
coefficients  in  the  representation  of  /  as  a  linear  combination  of  the  functions  P™,  P™+1, 
. . . ,  Pnh P^;  the  reconstruction  algorithm  computes  the  sample  values  of  /  from  the 
coefficients  in  the  representation  of  /  as  a  linear  combination  of  the  functions  Pm,  Pm+lT 
. . . ,  Pjv-i,  Pn-  The  CPU  times  are  in  seconds.  The  errors  are  relative  root-mean-square 
errors. 

The  columns  labeled  “‘fast’  transformation”  list  the  times  taken  by  the  algorithm  to 
evaluate  one  sum  of  the  form  (133),  (134),  (135),  or  (136). 

The  columns  labeled  “third  of  applying  an  N/2  x  N/2  matrix”  list  the  times,  divided 
by  3,  taken  to  apply  an  N/2  x  N/2  matrix  once  to  a  vector  of  length  N/2.  These  times 
scale  quadratically  with  N ;  the  figures  in  parentheses  are  estimates  used  when  the  memory 
required  by  direct  calculations  would  be  excessive. 

Remark  6.1  The  columns  labeled  “third  of  applying  an  N/2  x  N/2  matrix”  give  some  indi¬ 
cation  of  how  much  time  current  implementations  of  the  decompositions  and  reconstructions 
take  to  run;  these  times  are  believed  to  be  reasonable  estimates  of  how  a  standard  package 
like  Spherepack  (see  [2])  would  perform,  using  the  “semi-naive”  algorithm  described  in  [10], 
for  example.  Clearly,  for  sufficiently  small  N,  the  “slow”  implementation  would  actually  run 
faster  than  the  “fast”  implementation. 

The  code  was  compiled  with  the  Lahey-Fujitsu  compiler,  with  optimization  flag  — o2, 
and  run  on  a  2.8  GHz  Intel  Pentium  Xeon  microprocessor  with  512  Kb  of  L2  cache. 


22 


No  effort  was  made  to  optimize  the  precomputations.  To  simplify  the  implementation, 
precomputations  that  take  0(N2)  operations  were  used,  even  though  the  techniques  de¬ 
scribed  in  [4]  lead  naturally  to  precomputations  that  would  take  only  0(N  log  A" )  opera¬ 
tions.  The  precomputations  were  run  to  yield  approximately  6  digits  of  accuracy,  running 
all  computations  (including  the  precomputations)  in  double  precision  arithmetic. 

Observation  6.2  Asymptotically,  the  algorithm  should  require  a  number  of  operations  pro¬ 
portional  to  iV(logiV)  log(l/e)  to  evaluate  one  sum  of  the  form  (133),  (134),  (135),  or  (136), 
where  e  ~  10-6.  The  times  in  the  columns  labeled  “‘fast’  transformation”  appear  to  be  con¬ 
sistent  with  this  estimate. 


7  Generalizations 

The  algorithm  of  this  paper  admits  a  number  of  generalizations  and  extensions.  The  list 
below  is  not  intended  to  be  exhaustive;  subjects  in  it  are  under  investigation,  and  will  be 
reported  at  a  later  date. 

1.  Different  discretizations  of  S2.  Throughout  this  paper,  we  have  assumed  that  the 
meridians  on  S2  are  discretized  in  an  equispaced  manner,  i.e.,  that  the  points  60,  9i,  . . . , 
On- i,  On  subdivide  the  interval  [0,  7r]  into  equal  subintervals.  This  limitation  is  easily  re¬ 
moved  via  techniques  described  in  [11],  [16],  and  [5];  however,  to  maintain  numerical  stability, 
the  nodes  have  to  satisfy  certain  quite  restrictive  criteria.  One  important  collection  of  nodes 
that  does  in  fact  lead  to  stable  algorithms  is  given  by  the  formula 

0k  =  cos~1(xk),  (137) 

where  x0,  aq,  . . . ,  xn-i,  %n  are  Gaussian  nodes  on  the  interval  [—1,  1]  (see,  for  example,  [2]). 

Furthermore,  the  nodes  (p0,  </q,  . . . ,  </?2jv-i,  T2 n  in  the  discretizations  of  the  parallels  on 
S2  do  not  have  to  be  equispaced,  provided  some  form  of  the  non-equispaced  Fast  Fourier 
Transform  is  used;  again,  numerical  stability  requires  that  p 0,  </q,  . . . ,  P2N-1,  T2 n  be  fairly 
close  to  being  equispaced. 

2.  Associated  Laguerre  functions.  Associated  Laguerre  functions  are  defined  on  the 
entire  the  half-line  [0,  00),  but  are  very  small  outside  of  a  finite  interval.  The  Sturm-Liouville 
problem  that  generates  the  associated  Laguerre  functions  becomes  an  eigenvector  problem  for 
the  sum  of  a  diagonal  matrix  and  a  semiseparable  matrix  when  discretized  using  associated 
Laguerre  functions  of  low  orders  (very  much  like  the  associated  Legendre  functions). 

3.  Prolate  Spheroidal  Wave  Functions.  The  Sturm-Liouville  problem  that  generates  the 
Prolate  Spheroidal  Wave  Functions  becomes  an  eigenvector  problem  for  a  tridiagonal  matrix 
when  discretized  using  Legendre  polynomials. 

4.  Associated  Prolate  Spheroidal  Wave  Functions.  The  Sturm-Liouville  problem  that 
generates  the  associated  Prolate  Spheroidal  Wave  Functions  becomes  an  eigenvector  problem 
for  the  sum  of  a  tridiagonal  matrix  and  a  semiseparable  matrix  when  discretized  using 
associated  Legendre  functions  of  low  orders. 
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Table  1:  Times  in  seconds  and  errors  for  m  =  — even  degrees,  reconstruction 


N 

m 

“fast”  trans¬ 
formation 

third  of  applying 
an  f  xf  matrix 

precomp¬ 

utation 

relative 
r.m.s.  error 

1024 

515 

.20E-02 

. 91E-03 

. 20E+01 

. 64E-07 

2048 

1027 

.52E-02 

.36E-02 

. 80E+01 

. 58E-07 

4096 

2051 

.  23E-01 

. 15E-01 

. 33E+02 

. 17E-06 

8192 

4099 

.  48E-01 

.58E-01 

. 16E+03 

. 39E-06 

16384 

8195 

.  11E+00 

. 23E+00 

. 69E+03 

.72E-06 

32768 

16387 

. 24E+00 

. 96E+00 

. 36E+04 

. 14E-05 

65536 

32771 

.  51E+00 

( . 37E+01) 

. 14E+05 

. 52E-06 

131072 

65539 

.  11E+01 

( . 15E+02) 

. 63E+05 

.70E-05 

Table  2:  Times  in 

seconds  and  errors  for  m  =  N~(i ,  even  degrees,  decomposition 

N 

m 

“fast”  trans¬ 
formation 

third  of  applying 
an  £x£  matrix 

precomp¬ 

utation 

relative 
r.m.s.  error 

1024 

515 

. 19E-02 

. 91E-03 

. 30E+01 

.73E-08 

2048 

1027 

. 49E-02 

.36E-02 

. 80E+01 

.70E-08 

4096 

2051 

.  22E-01 

. 15E-01 

. 32E+02 

.76E-08 

8192 

4099 

.  47E-01 

. 58E-01 

. 16E+03 

. 11E-07 

16384 

8195 

.  11E+00 

. 23E+00 

. 68E+03 

. 13E-07 

32768 

16387 

. 23E+00 

. 96E+00 

. 35E+04 

. 15E-07 

65536 

32771 

. 50E+00 

( . 37E+01) 

. 14E+05 

.20E-07 

131072 

65539 

. 12E+01 

( . 15E+02) 

. 62E+05 

.35E-07 
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Table  3:  Times  in  seconds  and  errors  for  m  =  even  degrees,  reconstruction 


N 

m 

“fast”  trans¬ 
formation 

third  of  applying 
an  f  xf  matrix 

precomp¬ 

utation 

relative 
r.m.s.  error 

1024 

510 

.20E-02 

. 91E-03 

. 20E+01 

.77E-07 

2048 

1022 

.53E-02 

.36E-02 

. 90E+01 

.72E-07 

4096 

2046 

.  23E-01 

. 15E-01 

. 38E+02 

. 65E-07 

8192 

4094 

.  49E-01 

.58E-01 

. 17E+03 

. 15E-06 

16384 

8190 

.  11E+00 

. 23E+00 

. 78E+03 

. 41E-06 

32768 

16382 

. 24E+00 

. 96E+00 

. 34E+04 

. 54E-06 

65536 

32766 

. 52E+00 

( . 37E+01) 

. 15E+05 

. 31E-06 

131072 

65534 

.  11E+01 

( . 15E+02) 

. 64E+05 

.41E-05 

Table  4:  Times  in 

seconds  and  errors  for  m  =  N  2  4 ,  even  degrees,  decomposition 

N 

m 

“fast”  trans¬ 
formation 

third  of  applying 
an  £x£  matrix 

precomp¬ 

utation 

relative 
r.m.s.  error 

1024 

510 

. 19E-02 

. 91E-03 

. 20E+01 

.37E-08 

2048 

1022 

.  51E-02 

.36E-02 

. 90E+01 

.79E-08 

4096 

2046 

.  22E-01 

. 15E-01 

. 38E+02 

.95E-08 

8192 

4094 

.  48E-01 

. 58E-01 

. 16E+03 

. 11E-07 

16384 

8190 

.  11E+00 

. 23E+00 

. 77E+03 

. 16E-07 

32768 

16382 

. 24E+00 

. 96E+00 

. 33E+04 

. 15E-07 

65536 

32766 

.  51E+00 

( . 37E+01) 

. 16E+05 

. 17E-07 

131072 

65534 

.  11E+01 

( . 15E+02) 

. 62E+05 

.21E-07 
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Table  5:  Times  in  seconds  and  errors  for  m  =  jV3232 ,  odd  degrees,  reconstruction 


N 

m 

“fast”  trans¬ 
formation 

third  of  applying 
an  f  xf  matrix 

precomp¬ 

utation 

relative 
r.m.s.  error 

1024 

31 

. 16E-02 

. 91E-03 

. 20E+01 

.86E-08 

2048 

63 

.43E-02 

. 36E-02 

. 70E+01 

.20E-07 

4096 

127 

.  20E-01 

.  15E-01 

. 28E+02 

. 16E-06 

8192 

255 

.  43E-01 

.58E-01 

. 11E+03 

.41E-06 

16384 

511 

.  94E-01 

. 23E+00 

. 56E+03 

.60E-06 

32768 

1023 

.  21E+00 

. 96E+00 

. 26E+04 

. 10E-05 

65536 

2047 

.  46E+00 

( . 37E+01) 

. 14E+05 

. 14E-05 

131072 

4095 

.  10E+01 

( . 15E+02) 

. 47E+05 

. 14E-05 

Table  6:  Times  in  seconds  and  errors  for  m  =  ,V,J232 ,  odd  degrees,  decomposition 


N 

m 

“fast”  trans¬ 
formation 

third  of  applying 
an  f  xf  matrix 

precomp¬ 

utation 

relative 
r.m.s.  error 

1024 

31 

.15E-02 

.9 IE-03 

. 20E+01 

.71E-08 

2048 

63 

. 41E-02 

. 36E-02 

. 70E+01 

.59E-08 

4096 

127 

. 20E-01 

.  15E-01 

. 27E+02 

.82E-08 

8192 

255 

.  42E-01 

.58E-01 

. 12E+03 

.76E-08 

16384 

511 

.92E-01 

. 23E+00 

. 55E+03 

. 12E-07 

32768 

1023 

.  21E+00 

. 96E+00 

. 25E+04 

. 13E-07 

65536 

2047 

.  45E+00 

( . 37E+01) 

. 11E+05 

. 11E-07 

131072 

4095 

.  98E+00 

( .  15E+02) 

. 46E+05 

.20E-07 
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Table  7:  Times  in  seconds  and  errors  for  m  =  odd  degrees,  reconstruction 


N 

m 

“fast”  trans¬ 

third  of  applying 

precomp¬ 

relative 

formation 

an  f  xf  matrix 

utation 

r.m.s.  error 

1024 

960 

.  21E-02 

. 91E-03 

. 30E+01 

.12E-07 

2048 

1920 

.55E-02 

.36E-02 

. 11E+02 

. 20E-07 

4096 

3840 

. 23E-01 

. 15E-01 

. 39E+02 

. 48E-07 

8192 

7680 

. 49E-01 

. 58E-01 

. 16E+03 

. 14E-06 

16384 

15360 

. 11E+00 

. 23E+00 

. 83E+03 

. 15E-06 

32768 

30720 

.  24E+00 

. 96E+00 

. 36E+04 

. 23E-06 

65536 

61440 

.  53E+00 

( . 37E+01) 

. 21E+05 

. 12E-05 

131072 

122880 

.  11E+01 

( . 15E+02) 

. 62E+05 

.15E-05 

Table  8: 

Times  in 

seconds  and  errors  for  m  =  1  ,  odd  degrees,  decomposition 

N 

m 

“fast”  trans¬ 

third  of  applying 

precomp¬ 

relative 

formation 

an  matrix 

utation 

r.m.s.  error 

1024 

960 

.20E-02 

. 91E-03 

. 30E+01 

.48E-08 

2048 

1920 

.52E-02 

.36E-02 

. 10E+02 

.74E-08 

4096 

3840 

.  23E-01 

. 15E-01 

. 38E+02 

.10E-07 

8192 

7680 

. 49E-01 

. 58E-01 

. 16E+03 

. 11E-07 

16384 

15360 

. 11E+00 

. 23E+00 

. 81E+03 

.13E-07 

32768 

30720 

.  24E+00 

. 96E+00 

. 35E+04 

.18E-07 

65536 

61440 

. 51E+00 

( . 37E+01) 

. 16E+05 

. 16E-07 

131072 

122880 

.  11E+01 

( . 15E+02) 

. 60E+05 

.27E-07 
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